Alpha diversity
23 May 2024
1 Setup
1.1 Libs
## Loading required package: carData## Loading required package: lattice## Warning: package 'lattice' was built under R version 4.3.1## Loading required package: plyr## Warning: package 'plyr' was built under R version 4.3.1## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 4.3.1##
## Attaching package: 'ggpubr'## The following object is masked from 'package:plyr':
##
## mutate## Warning: package 'cowplot' was built under R version 4.3.1##
## Attaching package: 'cowplot'## The following object is masked from 'package:ggpubr':
##
## get_legend## Warning: package 'lme4' was built under R version 4.3.1## Loading required package: Matrix## Warning: package 'Matrix' was built under R version 4.3.1## Warning in check_dep_version(): ABI version mismatch:
## lme4 was built with Matrix ABI version 1
## Current Matrix ABI version is 0
## Please re-install lme4 from source or restore original 'Matrix' package## Warning: package 'glmmTMB' was built under R version 4.3.1## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.6.0
## Current Matrix version is 1.6.1.1
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package## Warning: package 'emmeans' was built under R version 4.3.1## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyser## Loading required package: stats4## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')##
## Attaching package: 'nlme'## The following object is masked from 'package:lme4':
##
## lmList##
## Attaching package: 'bestNormalize'## The following object is masked from 'package:MASS':
##
## boxcox1.2 Data
setwd("~/nicolagk@hawaii.edu - Google Drive/My Drive/Mosquito_business/Mosquito_microbes/04.microbiome_analysis/04a.diversity")
ps.clean <- readRDS("../../02.process_asvs/ps.clean.trim.less.rds")
ps.clean ## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 161 taxa and 458 samples ]
## sample_data() Sample Data: [ 458 samples by 24 sample variables ]
## tax_table() Taxonomy Table: [ 161 taxa by 11 taxonomic ranks ]2 Alpha diversity calculations (first run through)
##option to run rarefied stuff:
#ps.clean <- ps.rare
df.all <- data.frame(estimate_richness(ps.clean, split=TRUE, measures=c("Shannon","InvSimpson","Observed")))
df.all$sample_name <- rownames(df.all)
df.all.div <- merge(df.all,samdf,by="sample_name") #add sample data
#shannon diversity divided by species richness
df.all.div$even <- df.all.div$Shannon/(log(df.all.div$Observed))
df.div <- df.all.div
#adding most recent sample counts
sums <- data.frame(sample_sums(ps.clean))
colnames(sums) <- c("lib_size_trim")
sums$sample_name <- row.names(sums)
df.div <- merge(df.div,sums,by="sample_name")
#write.csv(df.div,"div.trim.csv")3 Analyzing alpha metrics
3.1 All data
Very strange/interesting pattern where the water types that shouldn’t have a lot of diversity have a lot of OTUs, but not a lot of reads…
ggplot(df.div,aes(x=type,y=lib_size_trim,color=infusion))+
geom_boxplot()+
geom_jitter(position=position_jitterdodge())
3.2 Plotz
#just experimental types
df.div.exp <- subset(df.div,type=="Microbial Water"|type=="A.albopictus")
df.div.exp$inf.temp <- paste0(df.div.exp$infusion,df.div.exp$temperature)
df.div.exp$type <- sub("A.albopictus","Mosquitoes",df.div.exp$type)
df.div.exp$type <- sub("Microbial Water","Meso. water",df.div.exp$type)
df.div.exp$type <- factor(df.div.exp$type,levels=c("Mosquitoes","Meso. water"))3.2.1 Richness
df.div.exp.se <- summarySE(df.div.exp,measurevar="Observed",groupvars=c("inf.temp","infusion","temperature","type"))
##dots & error bars
gg.rich <- ggplot(df.div.exp,aes(x=infusion,y=Observed,fill=inf.temp,shape=inf.temp))+
geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
#geom_boxplot(outlier.shape=NA,alpha=0.5)+
ylim(4,34)+
facet_wrap(~type)+
theme_cowplot()+
ylab("ASV richness")+
xlab("Infusion")+
scale_x_discrete(labels=c("OL","SG","PW"))+
scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))
gg.rich
##water by time
df.div.water <- subset(df.div.exp,type=="Meso. water")
gg.mw.rich.time <- ggplot(df.div.water,aes(x=infusion,y=Observed,fill=inf.temp))+
geom_boxplot()+
#geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
#geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
#geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
#geom_boxplot(outlier.shape=NA,alpha=0.5)+
facet_wrap(~day)+
theme_cowplot()+
ylab("ASV richness")+
xlab("Infusion")+
scale_x_discrete(labels=c("OL","SG","PW"))+
scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
ggtitle("Mesocosm water")
gg.mw.rich.time
3.2.2 Simpson’s
df.div.exp.se <- summarySE(df.div.exp,measurevar="InvSimpson",groupvars=c("inf.temp","infusion","temperature","type"))
##dots & error bars
gg.simp <- ggplot(df.div.exp,aes(x=infusion,y=InvSimpson,fill=inf.temp,shape=inf.temp))+
geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
geom_errorbar(data=df.div.exp.se,aes(ymax=InvSimpson+se,ymin=InvSimpson-se),width=0.2,color="black",position=position_dodge(width=0.6))+
geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
#geom_boxplot(outlier.shape=NA,alpha=0.5)+
facet_wrap(~type)+
theme_cowplot()+
ylab("Simpson's (inv.)")+
xlab("Infusion")+
scale_x_discrete(labels=c("OL","SG","PW"))+
scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
scale_shape_manual(values=c(22,24,22,24,22,24),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))
gg.simp
gg.mw.simp.time <- ggplot(df.div.water,aes(x=infusion,y=InvSimpson,fill=inf.temp))+
geom_boxplot()+
#geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
#geom_errorbar(data=df.div.exp.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
#geom_point(data=df.div.exp.se,size=2.5,position=position_dodge(width=0.6))+
#geom_boxplot(outlier.shape=NA,alpha=0.5)+
facet_wrap(~day)+
theme_cowplot()+
ylab("Simpson's index (inv.)")+
xlab("Infusion")+
scale_x_discrete(labels=c("OL","SG","PW"))+
scale_fill_manual(values=c("#B5D69F","#7DBA54","#C488A2","#992559","#B6F6EE","#60eede"),name="Infusion, temp.",labels=c("OL_cool","OL_warm","SG_cool","SG_warm","PW_cool","PW_warm"))+
ggtitle("Mesocosm water")
3.3 Statz
3.3.1 Richness
3.3.1.1 Water richness stats
df.div.mw <- subset(df.div.exp,type=="Meso. water")
df.div.mw$day <- as.factor(df.div.mw$day)
str(df.div.mw)## 'data.frame': 211 obs. of 30 variables:
## $ sample_name : chr "M_D13_OL1_1_S216_L001" "M_D13_OL1_2_S217_L001" "M_D13_OL1_3_S218_L001" "M_D13_OL1_4_S219_L001" ...
## $ Observed : num 25 28 22 28 31 31 28 29 27 21 ...
## $ Shannon : num 1.63 1.53 1.58 1.53 1.75 ...
## $ InvSimpson : num 2.67 2.4 2.44 2.32 2.93 ...
## $ orgname : chr "M_D13_OL1.1" "M_D13_OL1.2" "M_D13_OL1.3" "M_D13_OL1.4" ...
## $ mesocosm : chr "OL1.1" "OL1.2" "OL1.3" "OL1.4" ...
## $ type : Factor w/ 2 levels "Mosquitoes","Meso. water": 2 2 2 2 2 2 2 2 2 2 ...
## $ stage : chr "Water" "Water" "Water" "Water" ...
## $ sex : chr "" "" "" "" ...
## $ date.collected : chr "8/27/19" "8/27/19" "8/27/19" "8/27/19" ...
## $ day : Factor w/ 3 levels "4","12","20": 2 2 2 2 2 2 2 2 2 2 ...
## $ abdomen.length.mm: num NA NA NA NA NA NA NA NA NA NA ...
## $ wing.length : num NA NA NA NA NA NA NA NA NA NA ...
## $ hindleg.length : num NA NA NA NA NA NA NA NA NA NA ...
## $ special : chr "" "" "" "" ...
## $ dispersal : chr "N" "N" "D" "D" ...
## $ temperature : chr "C" "C" "C" "C" ...
## $ infusion : chr "OL" "OL" "OL" "OL" ...
## $ raw_miseq_reads : int 30391 39927 18829 58206 54860 39963 37857 45572 36841 27161 ...
## $ X16scq : num NA NA NA NA NA NA NA NA NA NA ...
## $ Wolb_ct1 : num NA NA NA NA NA NA NA NA NA NA ...
## $ Wolb_ct2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ Act_ct : num NA NA NA NA NA NA NA NA NA NA ...
## $ notes : chr "" "" "" "" ...
## $ newday : chr "early" "early" "early" "early" ...
## $ lib_size : num 25797 34139 15861 49298 44559 ...
## $ is.neg : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ even : num 0.506 0.458 0.511 0.46 0.51 ...
## $ lib_size_trim : num 25797 34134 15861 49290 44524 ...
## $ inf.temp : chr "OLC" "OLC" "OLC" "OLC" ...##
## Shapiro-Wilk normality test
##
## data: df.div.mw$Observed
## W = 0.9775, p-value = 0.001855#wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="poisson",data=df.div.mw)
#wrich2<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),family="compois",data=df.div.mw)
wrich3<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
#AICtab(wrich1,wrich2,wrich3)
##full model
wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
##no dispersal
wrich1nd<-glmmTMB(Observed~day+temperature+infusion+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus infusion
wrich1ni<-glmmTMB(Observed~day+temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus temperature
wrich1nt<-glmmTMB(Observed~day+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
##minus time
wrich1nnd<-glmmTMB(Observed~temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mw)
anova(wrich1,wrich1nd) #0.3931## Data: df.div.mw
## Models:
## wrich1nd: Observed ~ day + temperature + infusion + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1nd: (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1nd 8 1035.4 1062.2 -509.69 1019.4
## wrich1 9 1036.7 1066.8 -509.33 1018.7 0.7293 1 0.3931## Data: df.div.mw
## Models:
## wrich1ni: Observed ~ day + temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1ni: (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1ni 7 1118.0 1141.5 -552.02 1104.0
## wrich1 9 1036.7 1066.8 -509.33 1018.7 85.373 2 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.div.mw
## Models:
## wrich1nt: Observed ~ day + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1nt: (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1nt 8 1035.3 1062.2 -509.67 1019.3
## wrich1 9 1036.7 1066.8 -509.33 1018.7 0.6732 1 0.4119## Data: df.div.mw
## Models:
## wrich1nnd: Observed ~ temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1nnd: (1 | mesocosm), zi=~0, disp=~1
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1nnd 7 1045.2 1068.7 -515.61 1031.2
## wrich1 9 1036.7 1066.8 -509.33 1018.7 12.559 2 0.001874 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## day 12.9404 2 0.001549 **
## temperature 0.6742 1 0.411577
## infusion 151.4172 2 < 2.2e-16 ***
## dispersal 0.7306 1 0.392704
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Response: Observed
# Response: Observed
# Chisq Df Pr(>Chisq)
# day 12.9404 2 0.001549 **
# temperature 0.6742 1 0.411577
# infusion 151.4172 2 < 2.2e-16 ***
# dispersal 0.7306 1 0.392704
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
AICtab(wrich1,wrich1nd,wrich1ni,wrich1nt,wrich1nnd) #no temp best... almost tied with no dispersal, then full## dAIC df
## wrich1nt 0.0 8
## wrich1nd 0.1 8
## wrich1 1.3 9
## wrich1nnd 9.9 7
## wrich1ni 82.7 7##interactions?
wrich1.int <- glmmTMB(Observed~day*infusion+temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int1 <- glmmTMB(Observed~day+infusion*temperature+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int2 <- glmmTMB(Observed~day+infusion+temperature*dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int3 <- glmmTMB(Observed~day*temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int4 <- glmmTMB(Observed~day*dispersal+temperature+infusion+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.int5 <- glmmTMB(Observed~day+dispersal*infusion+temperature+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
anova(wrich1,wrich1.int) #2.86e-10***## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## wrich1.int: Observed ~ day * infusion + temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1036.66 1066.8 -509.33 1018.66
## wrich1.int 13 994.79 1038.4 -484.40 968.79 49.863 4 3.857e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## wrich1.int1: Observed ~ day + infusion * temperature + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int1: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1036.7 1066.8 -509.33 1018.7
## wrich1.int1 11 1037.2 1074.1 -507.62 1015.2 3.4154 2 0.1813## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## wrich1.int2: Observed ~ day + infusion + temperature * dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int2: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1036.7 1066.8 -509.33 1018.7
## wrich1.int2 10 1037.1 1070.7 -508.57 1017.1 1.5253 1 0.2168## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## wrich1.int3: Observed ~ day * temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int3: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1036.7 1066.8 -509.33 1018.7
## wrich1.int3 11 1039.0 1075.9 -508.52 1017.0 1.6226 2 0.4443## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## wrich1.int4: Observed ~ day * dispersal + temperature + infusion + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int4: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1036.7 1066.8 -509.33 1018.7
## wrich1.int4 11 1039.3 1076.2 -508.66 1017.3 1.3341 2 0.5132## Data: df.div.mw
## Models:
## wrich1: Observed ~ day + temperature + infusion + dispersal + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1: (1 | mesocosm), zi=~0, disp=~1
## wrich1.int5: Observed ~ day + dispersal * infusion + temperature + offset(log(lib_size_trim)) + , zi=~0, disp=~1
## wrich1.int5: (1 | mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1 9 1036.7 1066.8 -509.33 1018.7
## wrich1.int5 11 1037.8 1074.7 -507.90 1015.8 2.8599 2 0.2393## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## day 17.4811 2 0.00016 ***
## infusion 167.8537 2 < 2.2e-16 ***
## temperature 0.7595 1 0.38348
## dispersal 0.7627 1 0.38248
## day:infusion 58.9948 4 4.717e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: Observed
# Chisq Df Pr(>Chisq)
# day 17.4811 2 0.00016 ***
# infusion 167.8537 2 < 2.2e-16 ***
# temperature 0.7595 1 0.38348
# dispersal 0.7627 1 0.38248
# day:infusion 58.9948 4 4.717e-12 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#tapply(df.div.mw$log(Observed), df.div.mw$infusion, mean)
#tapply(df.div.mw$log(Observed), df.div.mw$day, mean)
wrich1.red <- glmmTMB(Observed~day+infusion+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
wrich1.red.int <- glmmTMB(Observed~day*infusion+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mw)
anova(wrich1.red,wrich1.red.int) #sig 4e-10***## Data: df.div.mw
## Models:
## wrich1.red: Observed ~ day + infusion + offset(log(lib_size_trim)) + (1 | , zi=~0, disp=~1
## wrich1.red: mesocosm), zi=~0, disp=~1
## wrich1.red.int: Observed ~ day * infusion + offset(log(lib_size_trim)) + (1 | , zi=~0, disp=~1
## wrich1.red.int: mesocosm), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wrich1.red 7 1034.08 1057.5 -510.04 1020.08
## wrich1.red.int 11 992.32 1029.2 -485.16 970.32 49.758 4 4.056e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Shapiro-Wilk normality test
##
## data: residuals(wrich1.red.int)
## W = 0.9921, p-value = 0.3128
#plotResiduals(wrich1.int)
#wrich1.red.int.resid <- simulateResiduals(fittedModel = wrich1.red.int, plot = T)
##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
wrich1.int.em <- emmeans(wrich1.int,~infusion*temperature)
multcomp::cld(wrich1.int.em)## infusion temperature emmean SE df lower.CL upper.CL .group
## SW C 19.1 0.351 198 18.5 19.8 1
## SW H 19.5 0.353 198 18.8 20.1 1
## SG C 21.2 0.355 198 20.5 21.9 2
## SG H 21.5 0.357 198 20.8 22.2 2
## OL C 24.7 0.351 198 24.0 25.4 3
## OL H 25.0 0.353 198 24.3 25.7 3
##
## Results are averaged over the levels of: day, dispersal
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.3.3.1.2 Water richness stats - rarefied
Same as above but without offset for library size. Not knitting this part, just checking out in R
df.div.mw <- subset(df.div.exp,type=="Meso. water")
df.div.mw$day <- as.factor(df.div.mw$day)
str(df.div.mw)
#hist(df.div.mw$Observed)
#shapiro.test(log(df.div.mw$Observed))
shapiro.test(df.div.mw$Observed)
#wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal)+(1|mesocosm),family="poisson",data=df.div.mw)
#wrich2<-glmmTMB(Observed~day+temperature+infusion+dispersal)+(1|mesocosm),family="compois",data=df.div.mw)
wrich3<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
#AICtab(wrich1,wrich2,wrich3)
##full model
wrich1<-glmmTMB(Observed~day+temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##no dispersal
wrich1nd<-glmmTMB(Observed~day+temperature+infusion+(1|mesocosm), data=df.div.mw)
##minus infusion
wrich1ni<-glmmTMB(Observed~day+temperature+dispersal+(1|mesocosm), data=df.div.mw)
##minus temperature
wrich1nt<-glmmTMB(Observed~day+infusion+dispersal+(1|mesocosm), data=df.div.mw)
##minus time
wrich1nnd<-glmmTMB(Observed~temperature+infusion+dispersal+(1|mesocosm), data=df.div.mw)
anova(wrich1,wrich1nd) #0.466
anova(wrich1,wrich1ni) #2.2e-16***
anova(wrich1,wrich1nt) #0.4789
anova(wrich1,wrich1nnd) #0.006693**
Anova(wrich1)
# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: Observed
# Chisq Df Pr(>Chisq)
# day 10.2547 2 0.005932 **
# temperature 0.5020 1 0.478604
# infusion 133.7700 2 < 2.2e-16 ***
# dispersal 0.5322 1 0.465695
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
AICtab(wrich1,wrich1nd,wrich1ni,wrich1nt,wrich1nnd) #no temp best... almost tied with no dispersal, then full
##interactions?
wrich1.int <- glmmTMB(Observed~day*infusion+temperature+dispersal+(1|mesocosm),data=df.div.mw)
wrich1.int1 <- glmmTMB(Observed~day+infusion*temperature+dispersal+(1|mesocosm),data=df.div.mw)
wrich1.int2 <- glmmTMB(Observed~day+infusion+temperature*dispersal+(1|mesocosm),data=df.div.mw)
wrich1.int3 <- glmmTMB(Observed~day*temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
wrich1.int4 <- glmmTMB(Observed~day*dispersal+temperature+infusion+(1|mesocosm),data=df.div.mw)
wrich1.int5 <- glmmTMB(Observed~day+dispersal*infusion+temperature+(1|mesocosm),data=df.div.mw)
anova(wrich1,wrich1.int) #4.35e-11***
anova(wrich1,wrich1.int1) #ns
anova(wrich1,wrich1.int2) #ns
anova(wrich1,wrich1.int3) #ns
anova(wrich1,wrich1.int4) #ns
anova(wrich1,wrich1.int5) #ns
Anova(wrich1.int)
# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: Observed
# Chisq Df Pr(>Chisq)
# day 13.9847 2 0.0009189 ***
# infusion 154.5860 2 < 2.2e-16 ***
# temperature 0.5795 1 0.4464988
# dispersal 0.5809 1 0.4459484
# day:infusion 64.8562 4 2.759e-13 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#tapply(df.div.mw$log(Observed), df.div.mw$infusion, mean)
#tapply(df.div.mw$log(Observed), df.div.mw$day, mean)
wrich1.red <- glmmTMB(Observed~day+infusion+(1|mesocosm),data=df.div.mw)
wrich1.red.int <- glmmTMB(Observed~day*infusion+(1|mesocosm),data=df.div.mw)
anova(wrich1.red,wrich1.red.int) #sig 4e-11**
##model checking
shapiro.test(residuals(wrich1.red.int)) #awesome
qqnorm(residuals(wrich1.red.int))
qqline(residuals(wrich1.red.int), col="red")
#plotResiduals(wrich1.int)
#wrich1.red.int.resid <- simulateResiduals(fittedModel = wrich1.red.int, plot = T)
##posthoc things
##just without dispersal because not plotting it
##also not plotting day but it should be in the model somewhere already
wrich1.int.em <- emmeans(wrich1.int,~infusion*temperature)
multcomp::cld(wrich1.int.em)3.3.1.3 Mosquito richness stats
##
## Shapiro-Wilk normality test
##
## data: df.div.mq$Observed
## W = 0.89464, p-value = 1.697e-10##different stat families
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
mrich2 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="compois",data=df.div.mq)
mrich3 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="genpois",data=df.div.mq)## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluationmrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="poisson",data=df.div.mq)
#mrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="nbinom1",data=df.div.mq)
#mrich5 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),family="nbinom2",data=df.div.mq)
AICtab(mrich1,mrich2,mrich3,mrich4) #gaussian again## dAIC df
## mrich1 0.0 10
## mrich2 47.1 10
## mrich4 48.3 9
## mrich3 48.3 10##full model
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##no dispersal
mrich1nd <- glmmTMB(Observed~newday+temperature+infusion+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus infusion
mrich1ni <- glmmTMB(Observed~newday+temperature+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus temperature
mrich1nt <- glmmTMB(Observed~newday+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
##minus time
mrich1nnd <- glmmTMB(Observed~temperature+infusion+dispersal+sex+offset(log(lib_size_trim))+(1|mesocosm), data=df.div.mq)
##minus sex
mrich1ns <- glmmTMB(Observed~newday+temperature+infusion+dispersal+offset(log(lib_size_trim))+(1|mesocosm),data=df.div.mq)
Anova(mrich1)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: Observed
## Chisq Df Pr(>Chisq)
## newday 1.8273 2 0.4011
## temperature 0.2244 1 0.6357
## infusion 0.4926 2 0.7817
## dispersal 1.4645 1 0.2262
## sex 0.3992 1 0.5275##
## Shapiro-Wilk normality test
##
## data: residuals(mrich1)
## W = 0.90866, p-value = 1.328e-09
## infusion temperature emmean SE df lower.CL upper.CL .group
## OL H 8.34 0.417 185 7.51 9.16 1
## SG H 8.49 0.447 185 7.61 9.37 1
## OL C 8.56 0.340 185 7.89 9.23 1
## SG C 8.72 0.490 185 7.75 9.68 1
## SW H 8.77 0.479 185 7.82 9.71 1
## SW C 8.99 0.564 185 7.88 10.10 1
##
## Results are averaged over the levels of: newday, dispersal, sex
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.3.3.1.4 Mosquito richness stats - rarefied
Same as above but without offset for library size
df.div.mq <- subset(df.div.exp,type=="Mosquitoes")
hist(df.div.mq$Observed)
shapiro.test(df.div.mq$Observed)
##different stat families
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),data=df.div.mq)
mrich2 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
mrich3 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="genpois",data=df.div.mq)
mrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="poisson",data=df.div.mq)
#mrich4 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="nbinom1",data=df.div.mq)
#mrich5 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="nbinom2",data=df.div.mq)
AICtab(mrich1,mrich2,mrich3,mrich4) #compois
##full model
mrich1 <- glmmTMB(Observed~newday+temperature+infusion+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
##no dispersal
mrich1nd <- glmmTMB(Observed~newday+temperature+infusion+sex+(1|mesocosm),family="compois",data=df.div.mq)
##minus infusion
mrich1ni <- glmmTMB(Observed~newday+temperature+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
##minus temperature
mrich1nt <- glmmTMB(Observed~newday+infusion+dispersal+sex+(1|mesocosm),family="compois",data=df.div.mq)
##minus time
mrich1nnd <- glmmTMB(Observed~temperature+infusion+dispersal+sex+(1|mesocosm),family="compois", data=df.div.mq)
##minus sex
mrich1ns <- glmmTMB(Observed~newday+temperature+infusion+dispersal+(1|mesocosm),family="compois",data=df.div.mq)
Anova(mrich1)
shapiro.test(residuals(mrich1))
qqnorm(residuals(mrich1))
qqline(residuals(mrich1), col="red")
mrich1.resid <- simulateResiduals(fittedModel = mrich1, plot = T)
##posthoc things
mrich1.em <- emmeans(mrich1,~infusion*temperature)
multcomp::cld(mrich1.em)
##no differences3.3.2 Simpson’s
3.3.2.1 Water Simpson’s stats
##
## Shapiro-Wilk normality test
##
## data: df.div.mw$InvSimpson
## W = 0.96756, p-value = 9.003e-05
##
## Shapiro-Wilk normality test
##
## data: log(df.div.mw$InvSimpson)
## W = 0.95673, p-value = 5.166e-06## Best Normalizing transformation with 211 Observations
## Estimated Normality Statistics (Pearson P / df, lower => more normal):
## - arcsinh(x): 1.5685
## - Box-Cox: 1.4779
## - Center+scale: 1.4303
## - Double Reversed Log_b(x+a): 1.5946
## - Exp(x): 3.4885
## - Log_b(x+a): 1.6482
## - orderNorm (ORQ): 1.2276
## - sqrt(x + a): 1.4479
## - Yeo-Johnson: 1.4679
## Estimation method: Out-of-sample via CV with 10 folds and 5 repeats
##
## Based off these, bestNormalize chose:
## orderNorm Transformation with 211 nonmissing obs and no ties
## - Original quantiles:
## 0% 25% 50% 75% 100%
## 1.359 1.947 2.669 3.225 5.130##
## Shapiro-Wilk normality test
##
## data: df.div.mw$simp.norm
## W = 0.99978, p-value = 1##just replaced InvSimpson with simp.norm
##full model
wsimp1 <- glmmTMB(simp.norm~day+temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##no dispersal
wsimp1nd <- glmmTMB(simp.norm~day+temperature+infusion+(1|mesocosm),data=df.div.mw)
##minus infusion
wsimp1ni <- glmmTMB(simp.norm~day+temperature+dispersal+(1|mesocosm),data=df.div.mw)
##minus temperature
wsimp1nt <- glmmTMB(simp.norm~day+infusion+dispersal+(1|mesocosm),data=df.div.mw)
##minus time
wsimp1nnd <- glmmTMB(simp.norm~temperature+infusion+dispersal+(1|mesocosm),data=df.div.mw)
Anova(wsimp1)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: simp.norm
## Chisq Df Pr(>Chisq)
## day 49.7901 2 1.542e-11 ***
## temperature 41.3001 1 1.306e-10 ***
## infusion 229.1055 2 < 2.2e-16 ***
## dispersal 0.0091 1 0.924
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: simp.norm
# Chisq Df Pr(>Chisq)
# day 49.7901 2 1.542e-11 ***
# temperature 41.3001 1 1.306e-10 ***
# infusion 229.1055 2 < 2.2e-16 ***
# dispersal 0.0091 1 0.924
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##~same results with normalized results
##same results with rarefied
AICtab(wsimp1,wsimp1nd,wsimp1ni,wsimp1nt,wsimp1nnd)## dAIC df
## wsimp1nd 0.0 8
## wsimp1 2.0 9
## wsimp1nt 32.6 8
## wsimp1nnd 40.7 7
## wsimp1ni 101.8 7##no dispersal good, then full model
##interaction?
wsimp1.int1 <- glmmTMB(simp.norm~day*infusion*temperature+(1|mesocosm), data=df.div.mw)
wsimp1.int2 <- glmmTMB(simp.norm~day*infusion+temperature+(1|mesocosm), data=df.div.mw) #best AIC
wsimp1.int3 <- glmmTMB(simp.norm~day+infusion*temperature+(1|mesocosm), data=df.div.mw)
wsimp1.int4 <- glmmTMB(simp.norm~day*temperature+infusion+(1|mesocosm), data=df.div.mw)
AICtab(wsimp1.int1,wsimp1.int2,wsimp1.int3,wsimp1.int4,wsimp1nd) #int better## dAIC df
## wsimp1.int2 0.0 12
## wsimp1.int1 5.6 20
## wsimp1.int3 40.5 10
## wsimp1nd 40.6 8
## wsimp1.int4 44.3 10## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: simp.norm
## Chisq Df Pr(>Chisq)
## day 70.733 2 4.370e-16 ***
## infusion 227.491 2 < 2.2e-16 ***
## temperature 39.606 1 3.107e-10 ***
## day:infusion 58.541 4 5.874e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Response: simp.norm
# Chisq Df Pr(>Chisq)
# day 70.733 2 4.370e-16 ***
# infusion 227.491 2 < 2.2e-16 ***
# temperature 39.606 1 3.107e-10 ***
# day:infusion 58.541 4 5.874e-12 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##all of the below was bad before normalizing
shapiro.test(residuals(wsimp1))##
## Shapiro-Wilk normality test
##
## data: residuals(wsimp1)
## W = 0.96148, p-value = 1.724e-05
## infusion temperature emmean SE df lower.CL upper.CL .group
## SG C -1.255 0.0934 202 -1.439 -1.0708 1
## SG H -0.656 0.0939 202 -0.841 -0.4712 2
## OL C -0.104 0.0927 202 -0.287 0.0786 3
## SW C 0.442 0.0927 202 0.259 0.6244 4
## OL H 0.495 0.0932 202 0.311 0.6783 4
## SW H 1.040 0.0932 202 0.857 1.2240 5
##
## Results are averaged over the levels of: day, dispersal
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same. # infusion temperature emmean SE df lower.CL upper.CL .group
# SG C -1.255 0.0934 202 -1.439 -1.0708 1
# SG H -0.656 0.0939 202 -0.841 -0.4712 2
# OL C -0.104 0.0927 202 -0.287 0.0786 3
# SW C 0.442 0.0927 202 0.259 0.6244 4
# OL H 0.495 0.0932 202 0.311 0.6783 4
# SW H 1.040 0.0932 202 0.857 1.2240 53.3.2.2 Mosquito Simpson’s
From manuscript pre-revisions: “In contrast, the Simpson’s index, which integrates evenness into a measure of alpha diversity, of the adult mosquito microbiome was lower in male mosquitoes relative to female mosquitoes (p<0.0001). This trend was associated with the increased dominance of wAlbB in males. The Simpson’s index of the adult mosquito microbiome did not vary significantly with the dispersal, aquatic chemistry, temperature, or time of emergence.”
##full model
msimp.all <- glmmTMB(simp.norm~newday+temperature+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without sex
msimp.nosex <- glmmTMB(simp.norm~day+temperature+infusion+dispersal+(1|mesocosm), data=df.div.mq)
##without dispersal
msimp.nodis <- glmmTMB(simp.norm~newday+temperature+infusion+sex+(1|mesocosm), data=df.div.mq)
##without infusion
msimp.noinf <- glmmTMB(simp.norm~newday+temperature+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without temperature
msimp.notem <- glmmTMB(simp.norm~newday+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)
##without time
msimp.notim <- glmmTMB(simp.norm~temperature+infusion+dispersal+sex+(1|mesocosm), data=df.div.mq)
Anova(msimp.all)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: simp.norm
## Chisq Df Pr(>Chisq)
## newday 2.3823 2 0.30387
## temperature 0.3062 1 0.58000
## infusion 2.7438 2 0.25363
## dispersal 2.8697 1 0.09026 .
## sex 67.4477 1 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Analysis of Deviance Table (Type II Wald chisquare tests)
#
# Response: simp.norm
# Chisq Df Pr(>Chisq)
# newday 2.2911 2 0.31805
# temperature 0.2987 1 0.58469
# infusion 2.7011 2 0.25909
# dispersal 2.8332 1 0.09234 .
# sex 67.5319 1 < 2e-16 ***
AICtab(msimp.all,msimp.nosex,msimp.nodis,msimp.noinf,msimp.notem,msimp.notim)## dAIC df
## msimp.notem 0.0 9
## msimp.notim 0.1 8
## msimp.noinf 0.4 8
## msimp.all 1.7 10
## msimp.nodis 2.4 9
## msimp.nosex 59.0 8#no time & no temp tied for best, then no infusion
##all of the below was bad before normalizing
shapiro.test(residuals(msimp.all))##
## Shapiro-Wilk normality test
##
## data: residuals(msimp.all)
## W = 0.98791, p-value = 0.0961
##posthoc things
msimp.all.em <- emmeans(msimp.all,~infusion*temperature)
multcomp::cld(msimp.all.em)## infusion temperature emmean SE df lower.CL upper.CL .group
## SG H -0.2061 0.152 185 -0.506 0.0939 1
## SG C -0.1150 0.167 185 -0.445 0.2150 1
## OL H -0.0782 0.143 185 -0.361 0.2045 1
## OL C 0.0129 0.115 185 -0.215 0.2406 1
## SW H 0.1484 0.162 185 -0.172 0.4682 1
## SW C 0.2394 0.193 185 -0.141 0.6197 1
##
## Results are averaged over the levels of: newday, dispersal, sex
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 6 estimates
## significance level used: alpha = 0.05
## NOTE: If two or more means share the same grouping symbol,
## then we cannot show them to be different.
## But we also did not show them to be the same.3.4 Infusion water things
df.div.inf <- subset(df.div,type=="Infusion water")
df.div.inf.se <- summarySE(df.div.inf,measurevar="Observed",groupvars=c("infusion"))
df.div.inf.se## infusion N Observed sd se ci
## 1 OL 3 45.666667 3.214550 1.8559215 7.985386
## 2 SG 3 23.666667 6.350853 3.6666667 15.776393
## 3 SW 3 4.333333 1.154701 0.6666667 2.868435# gg.iw.rich <- ggplot(df.div.inf,aes(x=infusion,y=Observed,fill=infusion))+
# geom_jitter(alpha=0.5,position=position_jitterdodge(jitter.width=0.2),color="gray")+
# geom_errorbar(data=df.div.inf.se,aes(ymax=Observed+se,ymin=Observed-se),width=0.2,color="black",position=position_dodge(width=0.6))+
# geom_point(data=df.div.inf.se,size=2.5,position=position_dodge(width=0.6))+
# #geom_boxplot(outlier.shape=NA,alpha=0.5)+
# theme_cowplot()+
# ggtitle("Infusion water")
# gg.iw.rich4 Diversity correlations
4.1 Diversity x day
Doesn’t look interesting
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
4.2 Diversity - mesocosm vs. mosquito
##
## Attaching package: 'reshape'## The following object is masked from 'package:Matrix':
##
## expand## The following object is masked from 'package:cowplot':
##
## stamp## The following objects are masked from 'package:plyr':
##
## rename, round_anymeso.corr <- merge(df.div.mq,df.div.mw,by="mesocosm")
ggplot(meso.corr,aes(x=InvSimpson.x,y=InvSimpson.y))+
geom_point()+
#facet_wrap(~day.y)+
geom_smooth()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'